

jDeltaEmissions<- function(delta.gen, plant.coef, summary=F, verbose=F){
  
  if(!'data.table'%in%class(delta.gen)) delta.gen = setDT(delta.gen)
  if('sample' %in% names(delta.gen)) delta.gen[,sample:=NULL]
  if(!'delta.NERC'%in%names(delta.gen)) setnames(delta.gen, old='pcac',new='delta.NERC')
  delta.gen[is.na(delta.gen)] = 0
  
  if(NROW(delta.gen)==0) return(NA)

    
  delta.gen.numeric = as.matrix(delta.gen[,grep('hr_', names(delta.gen)), with=F])/1000  #--> Solar data is in kW. Plant regressions have load in MW. Convert

  dimnames(delta.gen.numeric) = list(as.character(delta.gen$zip), grep('hr_', names(delta.gen), value=T))
  
  use.INT = unique(delta.gen[,'INT',with=F])
  use.NERC = unique(as.character(delta.gen[,delta.NERC]))
    if(NROW(use.INT)>1|NROW(use.NERC)>1) stop('Only one NERC/INT at a time!')
  use.ZIP = delta.gen[,'zip', with=F]
  
  if(use.INT=='TRE' ) cond.use = unique(plant.coef[grepl('_T',PCA),ORISPL])
  if(use.INT=='WEST') cond.use = unique(plant.coef[grepl('_W',PCA),ORISPL])
  if(use.INT=='EAST') cond.use = unique(plant.coef[grepl('_E',PCA),ORISPL])
  cond.use = as.character(cond.use)
  
  input = split(plant.coef[ORISPL%in%cond.use,], by = 'ORISPL', drop=T)
  
  res = lapply(input, function(zz, uN = use.NERC, dgn = delta.gen.numeric){ 
    zz = zz[grep(uN, zz$PCA),, drop=T][order(mo, hour, PCA),]
    if(NROW(zz)==0) return(NULL)
    pollcols = grep('_MASS', colnames(zz), value=T)
    ORISPL = unique(zz$ORISPL)
    stopifnot(NROW(ORISPL)==1)
    
    if(NROW(zz)!=288){ # fill in missing months with zeros
      eg = data.table(expand.grid(ORISPL = ORISPL, mo = 1:12, hour=min(zz$hour):max(zz$hour), PCA=zz$PCA[1]))
      zz = merge(zz, eg, by=c('ORISPL','hour','mo','PCA'), all=T)[order(mo, hour),][,(pollcols):=lapply(.SD, function(x) {x[is.na(x)] = 0   # 11-12-2019: DAMNIT! had: "function(x) is.na(x)=0" which returned all 0's for anything with an NA in it!
                                                                                                                          x}),.SD=pollcols]
    }
    
    zz = rbind(zz[rep(which(zz$mo==1),times=31),],
               zz[rep(which(zz$mo==2),times=28),],
               zz[rep(which(zz$mo==3),times=31),],
               zz[rep(which(zz$mo==4),times=30),],
               zz[rep(which(zz$mo==5),times=31),],
               zz[rep(which(zz$mo==6),times=30),],
               zz[rep(which(zz$mo==7),times=31),],
               zz[rep(which(zz$mo==8),times=31),],
               zz[rep(which(zz$mo==9),times=30),],
               zz[rep(which(zz$mo==10),times=31),],
               zz[rep(which(zz$mo==11),times=30),],
               zz[rep(which(zz$mo==12),times=31),])
    if(NROW(zz)!=8760) {
                        if(verbose==T) print(paste('Warning, ORISPL',zz$ORISPL_CODE[1],'does not have 8760 coefficients. Zeroing out', sep=" "))
                        dgn = rep(0, NROW(zz)) }
    
    tres = dgn%*%as.matrix(zz[,pollcols, with=F])
               
    data.frame(delta.zip = dimnames(delta.gen.numeric)[[1]],
               delta.ORISPL_CODE = ORISPL,
               delta.mWh = apply(delta.gen.numeric, MAR=1, sum, na.rm=T),
               delta.NERC = uN, 
               delta.CO2 = tres[,'CO2_MASS'], #--> CEM is in short tons, AP3 is in short tons. 
               delta.SO2 = tres[,'SO2_MASS']/2000, #--> CEM is in pounds, AP3 is in short tons: convert
               delta.NOX = tres[,'NOX_MASS']/2000,
               delta.PM2 = tres[,'PM2_MASS'], #--> input is in tons. All outputs are in tons! Except GLOAD
               delta.GLO = tres[,'GLOAD_MASS'], stringsAsFactors=F) 
  }) #--> end res lapply
  
  res = do.call(rbind, res)
  if(summary==T) return(aggregate(cbind(delta.CO2, delta.SO2, delta.NOX, delta.PM2, delta.GLO) ~ delta.zip + delta.SolarCap + delta.mWh + delta.NERC, data=res, FUN='sum', na.rm=T))
  if(summary==F) return(as.data.table(res))
}





jSummaryRes<-function(res, AP = AP3, cw = unique(Ot@data[,c('ORISPL_CODE','fips')]), detail=NULL){
  if(is.null(detail)) stop('Select detail=1,2,3 where 1 is lowest, 3 is highest')
  stopifnot('data.table'%in%class(res) & class(AP)=='array')
  res = res[,delta.ORISPL_CODE:=paste0('ORISPL_',delta.ORISPL_CODE)][order(res$delta.zip, res$delta.ORISPL_CODE),][which(res$delta.ORISPL_CODE%in%cw$ORISPL_CODE),]
  xres = cbind(res, cw[match(res$delta.ORISPL_CODE, cw$ORISPL_CODE),])
  xres = xres[,c('fips',grep(names(xres), pattern='CO2$|SO2$|NOX$|PM2$|GLO$', value=T)), with=F]
  
  #some xres's are NA's. When I dot-multiply, those NA's are dropped, which loops the non-NA's over again.
  # solve by setting NA's to 0
  stopifnot(sum(is.na(xres$fips))==0)
  xres[is.na(xres)] = 0
  
  
  ## Updated one FIPS: 12025 became 12086
  ## FIPS 26519 (in "res") does not exist. 26159 (MIchigan)
  xres[fips=='26519',fips:='26159']
  xres[fips=='12086',fips:='12025']
  xres = xres[, zip:=res$delta.zip]
  xres = xres[,j = lapply(.SD, sum), .SDcols = grep('delta', names(xres)), by=.(zip, fips)]
  # 3D array with rows equal to the FIPS in xres, columns are all fips, stack is pollutant
  xAP = AP3[xres$fips,,]
  dimnames(xAP)[[1]] = paste0(xres$zip, '--', xres$fips)
  
  # 3D array with rows equal to the FIPS in xres (plant locations), single column is the emissions in that fips (from panel in orig. zip), stack is pollutant 
  # need to rectify the different pollutants in AP3 vs. AP2
  #      dimnames(xAP)[[3]] # should match the order of the 3rd dim in xAP (which is the damages matrices).
  #      [1] "NOx"  "SO2"  "PM25" "CO2" 
  xx         = array(c(rep(xres$delta.NOX, dim(xAP)[2]),
                       rep(xres$delta.SO2, dim(xAP)[2]),
                       rep(xres$delta.PM2, dim(xAP)[2]),
                       rep(xres$delta.CO2, dim(xAP)[2])), 
                     dim=c(dim(xAP)[1],dim(xAP)[2],dim(xAP)[3]),
                     dimnames = dimnames(xAP)) * xAP
  if(detail==1){
    xx.sum = as.data.table(apply(xx, MARGIN=1, sum, na.rm=F))[,c('zip','fips'):=tstrsplit(dimnames(xx)[[1]], '--', keep=1:2)][,j = list(total_damages = sum(V1)), by=zip]
return(xx.sum)}
  
  if(detail==2){  # added this back 10/25/2019
    xx.sum = as.data.table(apply(xx, MARGIN=c(1,3), sum, na.rm=F))[,c('zip','fips'):=tstrsplit(dimnames(xx)[[1]], '--', keep=1:2)]
    xx.sum = xx.sum[,lapply(.SD, function(xx) sum(xx, na.rm=F)), by=list(zip), .SDcols = c('NOx','SO2','PM25','CO2')]
return(xx.sum)}
  if(detail==3){
    xx.sum = apply(xx, MARGIN=c(3), function(pp){
      as.data.table(pp)[,c('zip','fips'):=tstrsplit(dimnames(pp)[[1]], '--', keep=1:2)][,lapply(.SD, function(qq) sum(qq, na.rm=F)), by=list(zip), .SDcols = setdiff(colnames(pp), c('zip','fips'))]
    })
    return(rbindlist(xx.sum, idcol='POLL'))
  }
    } #--> returns a df with one row (named after zip) with total damages (if detail=1), POLL damages if detail=2, or 4 rows, 1 for each POLL + 1 col for each FIPS if detail=3




jSummaryResONELINE<-function(res, AP = AP3, cw = unique(Ot@data[,c('ORISPL_CODE','fips')]), detail=F){
  stopifnot(class(res)%in%c('data.table','data.frame') & class(AP)=='array')
  
  res = res[which(res$delta.ORISPL_CODE%in%cw$ORISPL_CODE),]
  # xres = merge(x=res, y=cw, by.x='delta.ORISPL_CODE', by.y='ORISPL_CODE', all.x=T, all.y=F, sort=F)
  xres = cbind(res, cw[match(res$delta.ORISPL_CODE, cw$ORISPL_CODE),])
  xres = xres[,c('fips',grep(names(xres), pattern='CO2$|SO2$|NOX$|PM2$|GLO$', value=T))]
  
  #some xres's are NA's. When I dot-multiply, those NA's are dropped, which loops the non-NA's over again.
  # solve by setting NA's to 0
  stopifnot(sum(is.na(xres$fips))==0)
  xres[is.na(xres)] = 0
  
  
  ## Updated one FIPS: 12025 became 12086
  ## FIPS 26519 (in "res") does not exist. 26159 (MIchigan)
  xres[which(xres$fips=='26519'),'fips'] = '26159'
  xres[which(xres$fips=='12086'),'fips'] = '12025'
  
  # 3D array with rows equal to the FIPS in xres, columns are all fips, stack is pollutant
  xAP = AP3[xres$fips,,]
  dimnames(xAP)[[1]] = row.names(res)
  
  # 3D array with rows equal to the FIPS in xres (plant locations), single column is the emissions in that fips (from panel in orig. zip), stack is pollutant 
  # need to rectify the different pollutants in AP3 vs. AP2
  #      dimnames(xAP)[[3]] # should match the order of the 3rd dim in xAP (which is the damages matrices).
  #      [1] "NOx"  "SO2"  "PM25" "CO2" 
  xres.array = array(c(rep(xres$delta.NOX, dim(xAP)[2]),
                       rep(xres$delta.SO2, dim(xAP)[2]),
                       rep(xres$delta.PM2, dim(xAP)[2]),
                       rep(xres$delta.CO2, dim(xAP)[2])), 
                     dim=c(dim(xAP)[1],dim(xAP)[2],dim(xAP)[3]),
                     dimnames = dimnames(xAP))
  
  
  # xres.array is emissions 
  # xAP is SAME SIZE, same names, is marginal damage per unit.
  # should be able to elementwise-multiply:
  xx = xres.array*xAP
  xx.sum = (apply(xx, MARGIN=1:2, sum, na.rm=F)) #--> MARGIN should *exclude* the dimension being summed over.
  if(detail==T) return(cbind(res, xx.sum[row.names(res),]))
  xx.sum2 = stats::aggregate(. ~ delta.zip + delta.SolarCap + delta.mWh + delta.NERC,
                             data= cbind(subset(res, select=-c(delta.ORISPL_CODE)), xx.sum[rownames(res),]),
                             sum)
  row.names(xx.sum2) = xx.sum2$delta.zip[1]
  return(xx.sum2)
} #--> returns a df with one row (named after zip) with total damages.




jPlotCoef<-function(OI, pollutant = 'CO2_MASS', lwd=0, coef.folder){
  OIn = OI
  if(class(OIn)!='numeric') OIn = as.numeric(sapply(strsplit(x= OI, split="_"), '[[', 2))
  if(class(OI)=='numeric') OI = paste0('ORISPL_',OI)
  
  ON = Ot$NERC[which(Ot$ORISPL==OI)]
  OS = Ot$UTLSRVNM[which(Ot$ORISPL==OI)]
  
  stopifnot(length(ON)>0)
  # cr = plant.coef[[OI]]
  cr = read.csv(file=file.path(coef.folder, paste0(OIn,'.csv')), stringsAsFactors=F)
  nms = cr$X
  
  cr$hour =       sapply(strsplit(x=nms, split=':'), '[[',1)
  cr$month =      sapply(strsplit(x=nms, split=':'), '[[',2)
  cr$delta.NERC = sapply(strsplit(x=nms, split='\\.'), '[[',2)
  
  cr$hour = as.numeric(sapply(strsplit(x=cr$hour, split='hour'),'[[',2))+1
  cr$month= as.numeric(sapply(strsplit(x=cr$month, split='mo'),'[[',2))
  
  cr$LineSize = as.numeric(cr$delta.NERC==ON)+lwd
  cr$pollutant = cr[,pollutant]
  
  return(ggplot(data=cr, aes(x=hour, y=pollutant, col=delta.NERC, size=LineSize)) + 
           geom_path() + facet_wrap(~month, nrow=2, ncol=6) +
           scale_size(range=c(1,1.75),guide=F) +
           labs(title=paste0(pollutant,' Marginal Response\n',OS,", ",ON, ", ", OI)))
}




